Estimates of multipolar coefficients to search 
for cosmic ray anisotropies with non-uniform 
or partial sky coverage. 
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Abstract 

We study the possibility to extract the multipolar moments of an underlying dis- 
tribution from a set of cosmic rays observed with non-uniform or even partial sky 
coverage. We show that if the degree is assumed to be upper bounded by L, each 
multipolar moment can be recovered whatever the coverage, but with a variance in- 
creasing exponentially with the bound L if the coverage is zero somewhere. Despite 
this limitation, we show the possibility to test predictions of a model without any 
assumption on L by building an estimate of the covariance matrix seen through the 
exposure function. 



1 Introduction 
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Anisotropy in the arrival directions of cosmic rays is a major observable to 
understand their origin. Magnetic fields bend their trajectories in such a way 
that transport of cosmic rays is mainly diffusive up to high energies: this 
makes their angular distribution isotropic. Nevertheless, above the so-called 
knee of cosmic rays up to the ankle, there are predictions for small but in- 
creasing anisotropies with energy, predictions which of course depend on the 
regular and the turbulent components of the assumed galactic magnetic field, 
as well as the assunaed distribution of so urces and composition of cosmic rays 
( iPtuskin et al. 19931 : ICandia et al. 20021 ). Further, at ultra-high energies, cos- 
mic ray arrival directions are expected to be less and less smeared out by galac- 
tic and extragalactic magnetic fields, leadin g to a possible ex t raction of infor 



mations about the pos i tion of the source s (llsola et al. 200ll: ISigl et al. 2003 



Armengaud et al. 20051 : iDolag et al. 2004 : iDe Marco et al. 20061 ). Hence, it is 
clear that any evidence for an anisotropy, or any limit on anisotropies in the 
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cosmic ray locations observed by experiments are among the most important 
constraints upon models. 



The multipole expansion up to a given order L is a powerful tool to study the 
structures standing out the noise down to an angular scale ^ vr/L, whatever 
the shape of the underlying celestial pattern. In practice, the number of sig- 
nificant coefficients is limited by the angular resolution of the detector and, 
in the other hand, by the available statistics of observation. However, ground 
based experiments cover a limited range in declination, so that it is impossible 
to apply off the shelf the formalism of multipole moments: anyone of the coeffi- 
cients may be modified in an unpredictable way by the unseen part of the sky. 
Methods have been developped to study the CMB with an incomplete coverag e 



(IGorski 19941 : IWright et al. 1994 iTeemark et al. 19961 : iMortlock et al. 20021 ) 



but here we are faced to a different problem: we cannot suppose a priori that 
the distribution of cosmic rays is described by a power spectrum, because we 
want to detect possible non-isotropic structures, a priori unknown. In other 
terms, the information carried by the aim cannot be reduced to the only knowl- 
edge of the Ci. 

One purpose of this paper is to study the possibility of estimating the mul- 
tipole moments of a distribution of points over a sphere in case of a non- 
uniform or even a partial coverage of the sky, together with the limitations of 
such an approa c h. The estimation of dipoles and quadrupoles was studied in 



(ISommers 200ll : lAubhn fc Parizot 20051 : iRoulet fc MoUerach 20051 1. Here, we 



use the moments of the observed distribution on a set of orthogonal functions: 
either the spherical harmonics themselves, or a set of functions tailored on 
the coverage function. With these two different methods, we show that the 
interference between the modes induced by the the non-uniformity or the hole 
of the coverage can be removed assuming a bounded expansion in the conju- 
gate space, allowing to recover the underlying multipole moments. However, 
in accordance with the simple intuition that it is impossible to describe the 
unseen part of the sky, we point out that the uncertainty on the recovered 
coefficients increases with the assumed bound L of the expansion. We show 
that the larger the hole in the coverage of the sky, the faster the increase of 
uncertainty with L. After some general considerations about the description 
of point processes on a sphere in Section 2, Sections 3 and 4 are dedicated to 
these methods whereas Section 5 illustrates them with some examples. 

Because of the incomplete knowledge of the distribution of cosmic ray sources, 
and the stochastic nature of the propagation through magnetic fields, the 
anisotropies we want to characterize are not reducible to explicit models: they 
may be interpreted as a particular realization of a random process. This means 
that some model predictions are better expressed as average values of the co- 
efficients, with their covariance matrix. This matrix is not necessarily diagonal 
to describe the physics we are interested in, contrary to the case of a power 
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spectrum. We show in Sect. 6 that under reasonable assumptions, an estimate 
can be performed with a partial sky coverage, evading the problem of setting 
a bound to the expansion. 



2 Generalities about point processes on a sphere 

The number of cosmic rays n{6, ip) observed as a function of = {9, ip) is a 
random process that we can modelize with the following quantity : 

1 ^ 

1=1 

where 5 is the Dirac function on the surface of the unit sphere, and the 
position of the i*'^ cosmic ray. This distribution follows a Poisson law with an 
averaged density that we will denote by ii{fl): 

Here, A is the density of the distribution of cosmic rays and lu is the exposure 
function of the experiment. The multipole coefficients of the function A(^, cp) 
defined on the unit sphere express its expansion in spherical harmonics: 

In this paper, we choose to normalize the spherical harmonics in such a 
way that / dQY£m{^)yi'm'{^) — ^T^^w^mm' ■ Together with the normalization 
/ dilA(Q) = 47r, our convention leads to ooo = 1 which is, in the context of 
this study, a natural system of units. For convenience, we will use hereafter 
the notation Y.i,m = E^o TL=-t 

With a uniform sky coverage, it is easy to obtain an unbiased evaluation of 
these coefficients from a sample of N points {6i, (pi) distributed independently 
according to the density A : 

1 ^ 

a^m = TrE^ri^i'^^O- 



If the distribution is roughly uniform (that is, \aim\ 1 for all {£, m) ^ (0, 0)), 
these estimators are quasi-optimal, weakly correlated and their variances are 
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Fig. 1. Relative exposure as a function of declination in equatorial coordinates of the 
Southern site of the Pierre Auger observatory. The detection efficiency is assumed 
to be saturated up to zenithal of 60° . 

close to otherwise the variances can be approximated from the quadratic 
moments: 



N 



[a in 



These properties are due to the orthogonahty of the spherical harmonics, and 
cannot been used directly if the coverage of the sphere is not uniform, that is, if 
the distribution actually observed is X{9, ip)u{9, (p), where a; is a non-uniform 
function eventually vanishing in some regions. 



However, if we suppose that the expansion of A in spherical harmonics is 
bounded to degree L (at least in good approximation), we are going to see 
that it is possible to recover - within limitations that we will discuss in details 
- the multipolar coefficients even in case of partial sky coverage. 
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Throughout the paper, we will consider by default an exposure function not 
covering t he whole sky in a realistic way since we use the function calculated by 
Sommers (jSommers 20011 ) describing the coverage of the sky of the Southern 
site of the Pierre Auger observatory as long as the acceptance of the detector 
is saturated until a local zenith angle ^^max- This function is shown on Fig.l 
with ^max = 60°, which guarantees in a realistic way this ideal function to be 
meaningful (lAuger Collaboration 20051 ). 



3 Estimate through the deconvolution of the exposure function 



3.1 The estimate 



In this section, we describe an estimate of the a^m coefficients based on the 
interpretation of the estimate 

1 ^ 

i=l 



in terms of a convolution between the underlying a^m coefficients of the density 
A(0, ip) and a kernel which depends on the uj{d, (p) function. In some extent, this 
approach is the eq uivalent of the MASTER one within the CMB framework 



( iHivon et al. 20021 ). except that we are interested here in building a linear 
estimate of the aim coefficients rather than a quadratic estimate of the Ci 
ones. As the cosmic rays are observed through the exposure function uj, the 
estimate h^rn is not an estimate of the multipolar coefficients of the density A, 
but an estimate of the multipolar coefficients of uj\. The a^rn coefficients are 
thus related to the ones through the following convolution 

I'm' 



The kernel K is entirely determined by the specific exposure function. Indeed, 
by using the completeness relation of the spherical harmonics, the elements of 
the kernel [K]^^?^' read 

4n 



This relation was refered to as the convolution theorem in ( iPeebles 19731 ). as 
this is the analog on the sphere of the convolution theorem for a Fourier's 
transform. Then, by using direct numerical results of K and for specific 
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exposure function u, the underlying a^* coefficients can be formally recovered 
through the following estimate 

e',m' 

3.2 Statistical properties of the estimate 

The observed points are sampled according to a Poissonian process on the 
sphere. Averaged over a large number of realisations of N events distributed 
according to fi{9,(f), it's elementar to compute the first and the second mo- 
ment oi n{fl) : 

(n(0))p = ^(0) 

{n{n)n{n'))p = ii{n)ii{n') + ii{n)d{n, n') 

where the subscript P stands for Poisson. The average of the estimate 
then reads 

{bem)p = (j dn n{n)Yr{n)\ 

Vtt ' p 

= y dQ i^{Q)Yp{Q) 

477 

i',m' 

leading to the following averaged a^m estimate 

li,mi h,m2 
— Clem- 

Thus, it is clear that we have built an unbiased estimate. Turning to the 
covariance, we get in the same way 

4lT 

4n 
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The only non vanishing term comes from the Poissonian part of the second 
moment oi n{fl): 

Using the fact that we are in practice looking for small deviation per respect 
with isotropy as emphasized in the introduction (ie: a^m/ooo <^ 1); this ex- 
pression can be simplified to: 

leading to : 

var(a,^) = [i^-i],T«oo. 



Let's remind that K being proportional to the number of events, the standard 
deviation of the reconstructed coefficients is hence proportional to 1/\/N as 
expected. In case of a non-uniform hut full coverage of the sky, the complete- 
ness relation of the spherical harmonics easily allows to give the following 
analytical expression of the operator : 



In case of partial coverage, the spherical harmonics are no longer orthogonal, 
in such a way that the coefficients of only satisfy the expression 



E mn^'Xr' = / dn^^Yr^Y^'in) 



where Ail is the non-zero region of ou, and 

[0]7r' = / dnYr{n)Y,f{n). 



It is then obvious, in this latter case, that K'^ is invertible only if L is fi- 
nite, and that the coefficients of strongly depend on the assumed bound 

L, leading to an indetermination of each coefficient as L is increasing. This 
indetermination is nothing else but the mathematical traduction that it's im- 
possible to know the distribution of cosmic rays in the uncovered region of the 
sky. 
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4 Estimate through dedicated orthogonal functions 



In this section, we describe another way, more intuitive, to recover the un- 
derlying aim coefficients by applying the Gram-Schmidt procedure to the 

uj{6,Lp)Yf^{9,Lp) with £ < L, which allows to build orthogonal functions from 
the coverage function. Then, by applying the formalism of moments to these 
functions; the are obtained with linear combinations of these moments. 

4-1 Applying the Gram- Schmidt procedure 
The scalar product being defined as 

{f\g)^j^Jne,^)g(e,^)dn, 

the normalized spherical harmonics may be written as 



where the are the associated Legendre functions supposed here to be 
normalized: 

1 } 

-jprixfdx^i. 



In practical computations we use the real functions Y^{6,ip) for m = 0, and 
{\/2PY^{6^ (/?) cos(m(/?), \^ Pf^{9, (p) sm{m(f)} for 1 < m < £. For convenience 
we keep the notations with the Y^^ hereafter. 

Let us suppose first that u; is a function of 9 only (for example, if the coverage 
is uniform in right ascension). Then uY^ and ujY^'' are orthogonal if m 7^ m', 
and the orthogonalisation may be performed separately for each value of m, 
combining the uoY^^ with m < £ < L. li represents the function / after 

normalization, we just need to set, for a given m: 



Q|I^|+2 =-^(^-P|m|+2 ~ {Q"m\\'^ ^1^1+2) Q"m\ ~ l^-P|m|+2) Q"m\+l) 
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p-1 

k=0 

Then the normahzed functions defined on the sphere by: 

{9, ip)^Qf {cos 9) e'""^ 

are orthogonal to each other, and the subset of with \m\ < i < L generates 
the same subspace as the ouYl" with \m\ < £ < L. We can express them through 
a set of coefficients C]^,: 

e'=m 

If uj depends on both 6 and ip, the same procedure can be appfied, but the 
orthogonal functions are mixtures of different values of m, and there is no 
canonical way to obtain them; anyway it is possible to build a basis preserving 
the subset generated by < £ < L whatever L. For simplicity, we do not 
develop such a formalism here. In particular, as only a small dependence on 
(/? is expected in the case we arc interested in, it is possible to weight the 
events to account for this variation of the exposure as a function of the right 
ascension, and hence, the formalism applied here can be applied off the shelf. 

To illustrate the method. Fig. 2 displays the shape of the ujYJ^ and the Z^, and 
the triangular matrix of transformation C^, for L = 15, m = (in logarithmic 
scale), in the case of the coverage function displayed in Fig.l. One can see 
that the off-diagonal terms (in absolute value) grow rapidly well above 1 with 
£ and dominate over the diagonal ones (the coefficients for other values of m 
have a quite simular pattern). This strong "mismatch" between the uYf"^ and 
the is suggested by the shapes of the functions and may be understood 
qualitatively in the following way: for large values of £—\m\, the 9 dependence 
of 17" has £ — \m\ oscillations over the full interval [0, tt], while Z'^ has the 
same number of oscillations over the covered interval; this makes difficult a 
matching of Z^^ to functions like a;17" which have less oscillations over this 
interval. 



4-2 Estimating the multipole coefficients 

Points being distributed according to the density X{9,(p) (to be evaluated), 
and detected with a probability u!{9) (supposed to be known), the observed 
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Fig. 2. Transformation between the wY^ and the with the coverage function of 
Fig.l. Top: shape (as a function of the declination) for £ < 4 (wY^ on left side, 
on right side). Bottom: coefficients for £ < 15; left: the C^£i, in logarithmic scale: 
the area (in the units of the axes) is In |C^^, |/10 (that is: a point represents 1, a unit 
square represents 2.2 x lO^j; the off-diagonal coefficients are in green if negative, 
in red if positive; right: the D^p (inverse matrix) in linear scale 1:1, with the same 
sign convention. 

points are distributed according to uX: this function may be expanded over 
the defined from uo as explained above: 



and an unbiased estimator of the aim is obtained from the points: 
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If A is quasi-uniform, ujX is almost proportional to Zq\ the coefficient aoo is 
largely dominant. Then these estimators are quasi-optimal; if uo is not constant, 
the a^mi foi' ^ given value of m, may be correlated. If A'^ is large, their covariance 
matrix is approximately given by quadratic moments: 

cov(a^^, aem) - T7 ^Ti^h Vi) Z^i^h Vi) - oiem Oiem- 



It is now easy to obtain estimators of the multipole coefficients of A at a given 
order L by substituting the expressions of the Z^^: 

that is: 

L L 

A ~ E ^^rn with aem = J2 ^t'm- 

e,m l'=l 

The 

a^jn with different values of m are not correlated, and the covariance 
matrix of the is given by: 

cov(a^i^,a^2^) = E Q™i C'f^fa cov(a^/^,a^"^). 



5 Illustrations 

To illustrate the statistical properties of the estimates, we show here some 
simple applications of the methods in case of exposure shown on Fig.l. For 
the sake of clarity, we will refer to as method 1 the method presented in section 
3, and to as method 2 the method presented in section 4. 

5. 1 Behaviour of variances with L 

For illustrations, we use here the method 1. In a first time, we restrict the 

bound L to 1, so that we are interested here in research of a dipolar compo- 
nent only. We show on Fig. 3 the reconstruction of the coefficient aio in the 
case of an indeed dipolar distribution, whose excess of events points towards 
equatorial North with a magnitude Oio = 0.1. The red histogram drawn show 
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Fig. 3. Reconstruction accuracy of the oiq coefficient in case of a aio = 0.1 dipolar 
pattern injected as a function of the assumed bound L=l (red), 2 (blue), 3 (green), 
in case of exposure function shown on Fig.l. Histograms are from Monte- Carlo, 
and superimposed curves are Gaussian with averages and standard deviations from 
analytical predictions. 

the occurence number of each reconstructed value of aio in case of = 10^ 
events generated by Monte-Carlo according to 



Over the histogram is plotted a Gaussian curve whose average and stan- 
dard deviation parameters are the ones determined in section 3.2. This curve 
matches the histogram, in such a way that the statistics previously deter- 
mined by calculation describe the properties of the estimators indeed. Let us 
note that under the assumption of a purely dipolar distribution (ie L=l) the 
reconstruction of the multipolar coefficients is obtained in a very reasonable 
way. 

Let us continue to illustrate the method by looking at the same multipolar 
coefficients, still in the case of a purely dipolar distribution, but by increasing 
the bound L to 2 and 3. Still on Fig. 3, the blue and the green histograms 
and Gaussian curves plot the same quantities than the red ones but for L = 2 
and L = 3 respectively, and illustrate the extremely fast degradation of the 
accuracy of the reconstruction of aio by more than a factor 2 for each additional 
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Fig. 4. Reconstruction accuracy of{a£m}£<2 coefficients as a function of the assumed 
bound L, in case of exposure function shown on Fig.l. Increasing L leads to an 
explosion of the uncertainty on aim ■ 



order. 



This tendency to the widening of the laws is largely confirmed when one looks 
at the reconstruction of any coefficients function of L. We show this 

property on Fig.4 for the {aem}e<2 set of coefficients, which illustrates clearly 
that it is increasingly difficult to give a meaning to the reconstructed values of 
the coefficients as soon as the maximum order of development is greater than 
3. 



5.2 Comparison of the two methods 



Two samples of points were simulated according to a slightly anisotropic dis- 
tribution (aio/aoo = 0.05, ai±i = 0, i.e. dipole moment along z), multiplied by 
the coverage function drawn in Fig. 1; the aim were estimated by both methods 
with the bound L going from 1 to 5. Fig. 5 shows that they give comparable 
results, and that the difference between them is generally smaller that the 
intrinsic difference between the samples (statistical fluctuations). Once again, 
one can see the divergence of the variances with increasing L. 
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Fig. 5. Comparison of reconstructed coefficients aim (i = 0,1) on two simulated 
dipolar samples (red, blue) with the two methods (solid circles: method 1 presented 
in section 3; open circles: method 2 presented in section 4 ), with a hound L from 1 
to 5. Top left: ooo; top right: oio; bottom: an and ai_i. 

5.3 Highly non-uniform coverarge of the whole sky 



A contrario, with a complete coverage (even highly non-uniform), the size of 
the variance is stabilized at large L. This is illustrated in Fig. 6, comparing a 
partial coverage (cf Fig. 1) to this coverage completed by a small fraction of the 
same function in the opposite hemisphere, in such a way that there is no fully 
unseen region. Even a relatively small relative exposure in the Northern part 
of the sky allows to recover the coefficients with almost the same precision 
as if the exposure was uniform on the whole sky. Note however that if the 
exposure in the opposite hemisphere tends to zero, even if the phenomenon 
of stabilization at large L remains, the variance at any i increases, tending 
towards a plateau determined roughly by where N' is in that case the 

total number of events which would be observed on the full sky through a 
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Fig. 6. Dependence on the bound L of the variances of estimated aim with a partial 
or a complete (but not uniform) coverage. Black circles: coverage function of Fig. 1; 
red squares: the same function plus 0.1 times the symmetric one w.r.t. the equatorial 
plane; blue triangles: the same + 0.2 times the symmetric one. The solid symbols 
correspond to a simulated dipolar distribution along z axis; the open ones correspond 
to a axis in the equatorial plane 

uniform window but with a low absolute coverage, in such a way that A^' is 
small. Of course, the larger the size of the relative exposure tending to is, 
the faster the increase of the variance towards this plateau occurs. 



5.4 Angular distribution in the covered region 



We have shown that using a large value of L in case of a partial coverage of 
the sky forbids to give to any a^m coefficient an interpretation of an individual 
multipolar moment. Nevertheless, one may wonder about the signification of 
the full set of coefficients {aim}- As a toy example, we generated a distribu- 
tion of points according to the exposure function of Fig.l times the function 
shown on Fig. 7 (top left) which is a combination of X2 ^"^^ ■ ^*^P 
right of Fig. 7, we show the reconstructed sky assuming L to be equals to 3, 
which illustrates that the reconstructed sky matches the injected one in the 
covered region even if the variance on each reconstructed multipolar coefficient 
is already large (as shown in preceding sub-sections) for L=3. Increasing the 
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Fig. 7. Top left: Toy injected sky (combination of ^^'^ 'Y^)' equatorial 

coordinates. Then recovered skies using different assumptions for L, and using the 
exposure function of Fig.l. Top right: Recovered sky assuming L=3. Bottom left: 
Recovered sky assuming L=5. Bottom right: Recovered sky assuming L=10. The 
unseen part of the sky is hidden. 

value of L to 5 (bottom left) or 10 (bottom right) do not change this property 
of the expansion, as only additional statistical fluctuations appear due to the 
finite number of points. On these plots, we hide the unseen part of the sky, 
where the reconstructed expansion is meaningless. 



5.5 Hypothesis test 

Any sky observed through an exposure function u can thus be described pre- 
cisely in the observed part of the sky by increasing L at a sufficient value. 
However the interpretation of each multipolar moment is problematic, be- 
cause it depends strongly on the cut L. We want now to build a statistical 
test to obtain a reasonable value of L from the data themselves. 

Starting from an hypothesis on L and the corresponding reconstructed {aem} 
coefficients, the likelihood function Cl built from the realization is 
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Fig. 8. Observed (red histograms) and expected (continuous blue line) distributions 
for A in case of a quadrupolar injected pattern. Top: L = 2 against the null hypothesis 
L=l. Bottom: L = 3 against the null hypothesis L=2. The vertical line indicates 
the value of A below which the null hypothesis is accepted with a threshold at 5%. 

For any particular realization, from this likelihood (which depends on L), 
we apply the method of the likelihood ratio to accept or to reject (within 
some chosen threshold) a null hypothesis Hq{Clq) with respect to another 
hypothesis Hi{Cl^) by computing 

n 



Asymptotically, for a sample obeying the hypothesis Hq, —2 In A is distributed 
according to a with a number of degrees of freedom equal to the number 
of extra parameters in the Hi hypothesis with respect to Hq. The value of A 
for any particular realization can thus be used to validate (or to reject) an 
assumption on L. 

As an example, let us assume that the cosmic rays follow a symmetrical 
quadrupolar distribution 1 + 0.1 sin^ 6* — 0.2 cos^ 6*, and let us use once again 
the exposure function shown on Fig.l. By restricting the reconstruction to a 
dipolar distribution, one then finds an artefact amplitude of about 5%. To test 
the relevance of the hypothesis of a purely dipolar sky, one can thus - starting 
from this sky - estimate whether it is necessary or not to increase the degree 



17 



of the expansion by calculating the ratio of the likelihood between the null 
hypothesis L=l and another hypothesis on L, L=2 for instance. To show the 
behavior of the test, we generated 1000 different realizations of the quadrupo- 
lar pattern with 100,000 points each, then we reconstructed the parameters 
of the expansion within the two hypotheses, and finally computed the ratio of 
likelihoods. In this case, the hypothesis L — 2 introduce 5 more parameters 
{a2m}, and the expected values of —2 In A are asymptotically distributed as a 
with 5 degrees of freedom. We plot the result on the top of Fig. 8: by choos- 
ing the threshold of the test to be 5% (vertical line at — 2 In A = 11.07), only 8 
reahzations over 1000 are accidentally accepted (red histogram). On the con- 
trary, repeating the same procedure to the L=2 null hypothesis with respect 
to the L=3 hypothesis, we show on the bottom of Fig. 8 that the obtained 
distribution perfectly matches the asymptotical expected one (a with 7 
degrees of freedom in that case). With the same partial coverage, a similar 
test on samples of 1000 points gives a poor discrimination between different 
hypotheses, and with only 100 points the test is completely irrelevant. 

This procedure may be used to define a "likely minimum value" Lj^in of L, 
and to prevent a wrong interpretation of multipolar coefficients obtained with 
a lower value, which are then biased (as the artefact dipole obtained above 
from a symmetric quadrupole). Of course, a given sample cannot provide by 
itself an absolute maximum for L, and in the presence of a hole the multipolar 
coefficients remain undefined without an external assumption; however let us 
point out that in many cases, the values of the coefficients at a given order 
have no intrinsic physical meaning if the distribution is of higher order. 



6 Testing model predictions 

Let us consider a distribution A(^, (/?) with coefficients on the Yf^; the 
observed distribution u!{9)X{9, ip) has coefficients aim on the Yf^, and the 
relation: 

oo 

dim = X] ^i^i '^^'"^ 

e'=e 

may be inverted, because for each value of m the matrix Cw is triangular, and 
the coefficients of the inverse relation may be computed exactly for any value 
of £ : 

oo 

e'=e 
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The values of the are displayed in Fig. 2 (right) with the same example as 
on the left plot, but in linear scale: contrary to the C™^, they remain below 1 
in absolute value, and practically negligible far away from the diagonal 

As a consequence, if a model gives predictions about the a^^, it will be possible, 
in some cases, to deduce predictions on the aim, which can be tested without 
any assumption on L. In that sense, the compatibility of a model may be 
checked with observations over an incomplete sky with a precision depending 
on the available statistics (but, of course, it can never prove that this model 
is the only possible one). 

If a model makes a deterministic prediction, comparing the aim to the pre- 
dicted values may be a convenient way to test this model up to a given order of 
multipolarity, that is, down to a given angular scale. The method is potentially 
more interesting if the predictions are probabilistic. As we emphasized it in 
the introduction, this is a relevant framework to describe high energy cosmic 
rays physics. Indeed, even in a situation with a well-defined and structured 
configuration of sources, propagation of cosmic rays unavoidably leads to a 
probabilistic nature of the obervable sky, that is to say, a probabilistic nature 
of the multipolar moments. Each class of models has intrinsically a natural 
variance encrypted in the a^m covariance matrix. Further, some models do 
not try to build a well-defined configuration of sources, but pick up randomly 
cosmic rays at sources according to some distributions, making even more im- 
possible to circumvent the characterization of a particular data set through a 
relevant statistical tool. 

Consequently, the discrimination of models through an exploratory search in 
a data set is potentially extremely powerful by looking for the distance of the 
full covariance matrix to the expected one. Most simple example is a model 
predicting random a^m following independent gaussian laws with variances aj. 
In that case, the covariance matrix for the a^m reads 

oo 

e'=m 



provided, of course, that this series converges: this is true if the series of o"| 
converges (as it should do for physical models), and if the D^^_^ have the 
behaviour suggested by Fig. 2. 



However, in practice, through the matrix inversion, the numerical divergence of 
the limits the expansion to L ~ 15 for this kind of coverage function; this is 
sufficient for most studies on sky anisotropies. 
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7 Conclusions 

To cope with a partial sky coverage, a formalism using the computation of 
moments on orthogonal functions was developed to recover the angular dis- 
tribution of the incident flux from a sample of observed points. If the mul- 
tipolar expansion is assumed to be upper bounded by £ < L, the coefficients 
aim may be estimated with a variance proportional to as usually, with a 
penalty factor increasing exponentially with L if there is a hole in the coverage 
(but stabilizing rapidly if the coverage is nowhere vanishing, even highly non- 
uniform). Two methods were tested, giving similar results, and practically the 
same variances. 

Statistical tests based on likelihood ratios may be built to check an hypothesis 
on the distribution, for example a given bound i < L. In any case, it is 
possible to express predictions of a model in terms of coefficients which can 
be computed without any assumption on L, and tested against the moments 
found with a sample of observed points. 

The methods presented in this paper may be applied any cosmic ray dataset, 
provided that the arrival directions and the coverage of the sky are known 
within a reasonable precision. 
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